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O : Abstract 

I The far-zone flux of energy contains hereditary (tail) contributions that depend on the entire past history of 

Q^' the source. Using the multipolar post-Minkowskian wave generation formalism, we propose and implement 

■^ ■ a semi-analytical method in the frequency domain to compute these contributions from the inspiral phase of 

a binary system of compact objects moving in quasi-elliptical orbits up to third post-Newtonian (3PN) order. 

The method explicitly uses the quasi-Keplerian representation of elliptical orbits at IPN order and exploits 

the doubly periodic nature of the motion to average the 3PN fluxes over the binary's orbit. Together with 

^1 the instantaneous (non-tail) contributions evaluated in a companion paper, it provides crucial inputs for the 

construction of ready-to-use templates for compact binaries moving on quasi-elliptic orbits, an interesting 

^ . class of sources for the ground based gravitational wave detectors such as LIGO and Virgo as well as space 

based detectors like LISA. 
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I. INTRODUCTION 

The gravitational- wave (GW) energy flux from a system of two point masses in elliptic mo- 
tion in the leading quadrupolar approximation (Newtonian order) was first obtained by Peters & 
Mathews lli |2|]. Using the first post-Newtonian (IPN) order quasi-Keplerian representation of 
the binary's orbit fS'], Blanchet & Schafer ^ computed the IPN corrections to the above result 
(confirming earlier work by Wagoner & Will 151]).' Using the generalized quasi-Keplerian rep- 
resentation of the 2PN motion la LZI, ISD, Gopakumar and Iyer ^ extended these results to 2PN 
order and computed the 'secular' evolution of orbital elements under 2PN gravitational radiation- 
reaction (4.5PN terms in the equations of motion). These constitute one of the basic inputs for 
gravitational wave phasing of binaries in quasi-eccentric orbits in the adiabatic approximation. 
All these works above relate to the instantaneous terms in the phasing of gravitational waves. 

The multipole moments describing GWs emitted by an isolated system do not evolve indepen- 
dently. They couple to each other and with themselves, giving rise to non-linear physical eff'ects. 
Consequently, starting at relative 1.5PN order, the above instantaneous terms in the flux must be 
supplemented by the contributions arising from these non-linear multipole interactions. The lead- 
ing multipole interaction is between the mass quadrupole moment /,y and the mass monopole M 
or Amowitt, Deser and Misner (ADM) mass. It is associated with the non-linear effect of tails at 
order 1 .5PN, and is physically due to the backscatter of linear waves from the space-time curva- 
ture generated by the total mass M. Tails imply a non-locality in time since they are described as 
integrals depending on the history of the source from the remote past to the current retarded time. 
They are thus appropriately referred to as hereditary contributions by Blanchet & Damour [llOl llin 



terms non-local in time depending on the dynamics of the system in its entire past llllll . The most 



detailed study of tails in this context 11121 11311 is based on the multipolar post-Minkowskian formal- 
ism lll4l ll5n. Up to 3PN order the hereditary terms comprise the dominant quadratic-order tails, 
the quadratic-order memory integral llllLll6Ul7Ul8U l9ll and the cubic-order tails. The latter cubic 



"monopole-monopole-quadrupole" interaction can be called "tails of tails" of GWs (see [|12l 11311 
for earlier references to the general topic of tails). In this paper we set up a general theoretical 
framework to compute the hereditary contributions for binaries moving in elliptical orbits and 
apply it to evaluate all the tail contributions contained in the 3PN accurate GW energy flux. 

For the instantaneous terms in the energy flux, explicit closed form analytical expressions can 
be given in terms of dynamical variables related to relative velocity and relative separation. Con- 
sequently, these expressions can be conveniently averaged in the time domain over an orbit using 
their quasi-Keplerian representation. For the hereditary contributions on the other hand one can 
only write down formal analytical expressions as integrals over the past. More explicit expressions 
in terms of the dynamical variables a priori require a model of the binary's orbital evolution in 
the past to implement the integration. In general, one can show [lllll that the past-influence of tails 
decreases with some kernel oc l/(t - t')^ where t is the current time and t' the integration time in 
the past. Thus the "remote-past" contribution to the tail integrals is negligible. More precisely. 



it was shown [|20ll that the contribution due to the past of the tail integral is (9(^rad In ^rad) where 
^lad = co/o/ is the adiabatic parameter associated with the binary's inspiral due to radiation re- 
action, which is of order 2.5PN. Consequently, the tail integrals may be evaluated using standard 
integrals for a fixed non-decaying circular orbit and errors due to inspiral by gravitation radiation 



reaction are at least 4PN order [|20ll . 



' As usual the «PN order refers to the post-Newtonian terms of order (u/c)^" where u denotes the typical binary's 
orbital velocity and c is the speed of hght. 



In the circular orbit case, with the above simplified model of binary inspiral one can work 
directly in the time domain. For instance, the hereditary terms in the flux were computed up to 



3.5PN [112|,11J] while those in the GW polarisations could be obtained up to 2.5PN [tl9i,l21l]. In 
the elliptic orbit case on the other hand the situation is more involved. Even after using the quasi- 
Keplerian parametrization, one cannot perform the integrals in the time domain (as for the circular 
orbit case), since the multipole moments have a more complicated dependence on time and the 
integrals are not analytically solvable in simple closed forms. By working in the Fourier domain 
Ref. r2Q1 computed the hereditary tail terms at 1.5PN for elliptical orbits using the lowest order 
Keplerian representation. 

In the present investigation we tackle the terms at orders 2.5PN and 3PN and we need to go 
beyond the (Newtonian) Keplerian representation of the orbit to a IPN quasi-Keplerian represen- 
tation. Here we encounter two further complications. Firstly, the IPN parametrization of the 
binary [SI] involves three kinds of eccentricities (e,-, e, and e^). More seriously, at IPN order the 
periastron precession eff'ect appears in the problem and one has to contend with two times scales: 
the orbital time scale and the periastron precession time scale. These new features are to be prop- 
erly accounted for in the calculations to extend the Fourier method of Ref. 12011 . This strategy 
has been proposed and used earlier in computing the instantaneous terms in the GW polarizations 



from binaries on elliptical orbits 11221. I23L 12411 . We shall adapt these features here to treat the more 
involved hereditary contribution to the total energy flux. 

Following Ref. EOll . we express all the multipole moments needed for the hereditary compu- 
tation at Newtonian order as discrete Fourier series in the mean anomaly of motion i. However, 
for the quadrupole moment /,y needed beyond the lowest Newtonian order, the "doubly periodic" 
nature of the motion needs to be crucially incorporated. The evaluation of the Fourier coefficients 
is done numerically based on a series of combinations of Bessel functions. All tail terms at 2.5PN 
and 3PN are computed to provide the "enhancement factors" (functions of eccentricity playing 
a role similar to the classic Peters & Mathews [jl|l enhancement factor) for binaries in elliptical 
orbits at the 2.5PN and 3PN orders. The present work extends results for hereditary contributions 
at 1.5PN [I20I1 for elliptical orbits to 2.5PN and 3PN orders. The 3PN hereditary contributions 



comprise the tail-of-tail terms and are also extensions of [Il2l ll3ll for circular orbits to the elliptical 



case.^ 



Combining the hereditary contributions computed in this paper with the instantaneous contri- 



butions computed in the companion paper 



2911 will yield the complete 3PN energy flux, general- 



izing the circular orbit results at 2.5PN Il30ll and 3PN I131LI32L 13311 to the elliptical orbit case. The 
final expressions represent GWs from a binary evolving adiabatically under gravitational radiation 
reaction, including precisely up to 3PN order the effects of eccentricity and periastron precession 
during epochs of inspiral when the orbital parameters are essentially constant over a few orbital 
revolutions. It thus represents the first input to go towards the full quasi-elliptical case, namely 
the evolution of the binary in an elliptical orbit under gravitational radiation reaction. 

Recently, Damour, Gopakumar & Iyer [|24ll proposed an analytic method based on an improved 
"method of variation of constants" to construct high accuracy templates for the GW signals from 
the inspiral phase of compact binaries moving in quasi-elliptical orbits. The three time scales, re- 
spectively related to orbital motion, orbital precession and radiation reaction, are handled without 
the usual approximation of assuming adiabaticity relative to the radiation reaction time scale. The 



explicit results of the above treatment Il24ll relate to "Newtonian" radiation reaction (2.5PN terms 



^ Recall that tails are not just mathematical curiosities in general relativity but facets that should showup in the GW 



signals of inspiraling compact binaries and be decoded by the detectors Virgo/LIGO and LISA 11251 1261 1271 12811 



in the equations of motion). It leads to post-adiabatic (fast) oscillations resulting in amplitude cor- 



rections at order 2.5PN beyond the secular terms. More recently this work has been extended Il34|] 
to IPN radiation reaction (3.5PN terms in the equations of motion).^ 

The paper is organized as follows: In Section |II] we review the solution of the equations of 
motion of compact binaries and discuss its important properties relevant for this present work. 
Section |III] provides the Fourier decomposition of multipole moments and its use in averaging 
the energy flux. Section ITVl provides the computations of all the tail contributions whose numer- 
ical implementation is elaborated in Section |Vl The complete 3PN contributions are exhibited in 
Section |Vl] together with relevant checks. The paper ends with an Appendix listing the Fourier 
coefficients of the required Newtonian moments in terms of the Bessel functions. 

II. SOLUTION OF THE EQUATIONS OF MOTION OF COMPACT BINARIES 
A. Doubly-periodic structure of the solution 



In this work and the next one 112911 . we shall often need to use the explicit solution for the 
motion of non-spinning compact binary systems in the post-Newtonian (PN) approximation. We 
review here the relevant material we need, which includes the general "doubly-periodic" structure 
of the PN solution, and the quasi-Keplerian representation of the IPN binary motion by means of 
different types of eccentricities. We closely follow the works [l3 U22l l35n. 

The equations of motion of a compact binary system up to the 3PN order admit, when neglect- 
ing the radiation reaction term at the 2.5PN order, ten first integrals of the motion corresponding to 



the conservation of energy, angular and linear momenta, andpositionof the center of mass Il36l l37l1. 
When restricted to the frame of the center of mass, the equations admit four first integrals associ- 
ated with the energy E and the angular momentum vector J, given at 3PN order by Eqs. (4.8)-(4.9) 



ofRef. 11381] 



The motion takes place in the plane orthogonal to J. Denoting by r = |x| the binary's orbital 
separation in that plane, and by v = Vi -V2 the relative velocity, we find that E and J are functions of 
r, r^, v^ and XX V (we are employing for definiteness the harmonic coordinate system of 113 811 "^1. and 
depend on the total mass m = mi + mj and reduced mass fi = mim2/m. We adopt polar coordinates 
r, (p in the orbital plane, and express E and the norm 7 = |J|, thanks to v^ = r^ + r^(p^, as some 
explicit functions of r, f- and (p. The latter functions can be inverted (by means of straightforward 
PN iteration) to give r^ and (p in terms of r and the constants of motion E and J. Hence, 

r^ = nr-E,Jl (2.1a) 

4> = g[r,E,J], (2.1b) 

where the functions 'R and Q denote certain polynomials in 1/r, the degree of which depends on 



the PN approximation in question (it is seventh degree for both 'R and at 3PN order 13911 ). The 
various coefficients of the powers of 1 /r are themselves polynomials in E and J, and also, of 
course, depend on m and the dimensionless reduced mass ratio v = /u/m. In the case of bounded 



3 



For circular orbits, secular evolution of the phase, computed in the adiabatic approximation, is known up to 3.5PN 



order 1321 133|]. 
^ All calculations in this paper will be done at the relative IPN order, and at that order there is no difference between 
the harmonic and ADM coordinates. 



elliptic -like motion, one can prove [12211 that the function H admits two real roots, r-p and ta such 
that r-p < rA, which admit some non-zero finite Newtonian limits when c ^> oo, and represent 
respectively the radii of the orbit's periastron and apastron. The other roots tend to zero when 
c ^> oo. 

We are considering a given binary's orbital configuration, fully specified by some given values 
of the integrals of motion E and J. We no longer indicate the dependence on E and 7 which is 
always implicit in what follows. The binary's orbital period, or time of return to the periastron, is 
obtained by integrating the radial motion as 



P = 2 f — ^. (2.2) 

J rp 



dr 



1 C'' 
K=- \ dr^^, (2.3) 



We introduce the fractional angle (i.e. the angle divided by 2n) of the advance of the periastron 
per orbital revolution, 

which is such that the precession of the periastron per period is given by Acp = 2n(K - 1). As 
K tends to one in the limit c ^ oo (as is easily checked from the Newtonian limit), it is often 
convenient to pose k = K - I, which will then entirely describe the relativistic precession. 
Let us define the mean anomaly i and the mean motion n by 

I = n{t-tp), (2.4a) 

In 
n = —. (2.4b) 

Here fp denotes the instant of passage to the periastron. For a given value of the mean anomaly £, 
the orbital separation r is obtained by inversion of the integral equation 



^ = n r -^==. (2.5) 



dr' 



This defines the function r{i) which is a periodic function in i with period In. The orbital phase (p 
is then obtained in terms of the mean anomaly £ by integrating the angular motion as 

1 r' 

= 0P + - de'QW)], (2.6) 

n Jo 

where 0p denotes the value of the phase at the instant tp. In the particular case of a circular orbit, 
r = const, the phase evolves linearly with time, ^ = Q [r] = co, where co is the orbital frequency of 
the circular orbit given by 

co = Kn = {\+k)n. (2.7) 

In the general case of a non-circular orbit it is convenient to keep the definition of a» = Kn (which 
will notably be very useful in the next work [1290) and to explicitly introduce the linearly growing 



part of the orbital phase ([2.6[) by writing it in the form 



(f> = (f)p + co{t-tp) + W(£) 
= (Pp + K£ + W{£). (2.8) 



Here W{£) denotes a certain function which is periodic in £ (hence, periodic in time with period 
P). According to (12.61) this function is given in terms of the mean anomaly £ by 

1 r*" 

W(£) = -l df[g[rin]-(4- (2.9) 

Finally, the decomposition (12.81) exhibits clearly the "doubly periodic" nature of the binary motion, 
in terms of the mean anomaly £ with period 2n, and in terms of the periastron advance K € with 
period 2n K. ^ It may be noted that in Refs. 11231 12411 the notation A is used; it corresponds to 
A = K€ and will also occasionally be used here. 



B. Quasi-Keplerian representation of the motion of compact binaries 



In the following we shall also use the explicit solution of the motion at IPN order, in the form 
due to Damour & Deruelle iJ] . The solution is given in parametric form in terms of the eccentric 
anomaly u. Then the radius r and mean anomaly £ are expressed as 

r = ar(l - fircosw), (2.10a) 

i = u- e,smu. (2.10b) 

The phase angle (p is given by (the additive constant cp^ is for convenience set equal to zero) 

(P = KV, (2.11) 

where the true anomaly V is defined by^ 



V = 2 arctan 



1 \l/2 

-. -] tan- 

1 - e^ / 2 



(2.12) 



In the above, K is the periastron advance given in general terms by Eq. (12.31) . and Qr is the semi- 
major axis of the orbit. Note that there are, in this parametrization at IPN order, three kinds of 
eccentricities Cr, e, and e^ (labelled after the coordinates r, t and 0). All these eccentricities differ 
from one another by IPN terms, while the advance of the periastron per orbital revolution appears 
also starting at the IPN order. Due to these features, this representation is referred to as the "quasi- 
Keplerian" (QK) parametrization for the IPN orbital motion of the binary. The periodic function 
W of Eq. (12.91 ) now reads 

W = K{V-i). (2.13) 

To close the above solution we need to know the explicit dependence of the orbital elements in 
terms of the IPN conserved energy E and angular momentum J in the center-of-mass frame (taken 



^ Recall, that though standard, the term "doubly periodic" may mislead a bit in that the motion in physical space is 
not periodic in general. The radial motion r{t) is periodic with period P while the angular motion (pit) is periodic 
[modulo In} with period Pjk where k - K - \. Only when the two period are commensurable, i.e. when k = l/N 
where A^ is a natural integer, is the motion periodic in physical space (with period NP). 

* We have denoted the true anomaly by V rather than by the symbol u of earUer papers to avoid conflict with the 
relative speed v. 



as usual per unit of the reduced mass fi). This is given in Ref. yD. Note that the semi-major axis 
Gr and mean motion n depend at IPN order only on the constant of energy through 



,,_2^j,4^_r\£j. 



n = 



IE 
Gm 



2 1 c 



1 I 15 v\E 



Posing h = J/{Gm), the IPN periastron precession simply reads^ 

3 



K= 1 + 



c^h?-' 



while the three different eccentricities are given by 
e,. = \\+2Eh^ 
l+lEh^ 



,15 5 \E -6 + v 



1/2 



^ ,17 7 \E 2-2v 

.15 v\E 6 

2 ^2/c2 c2;z2 



1/2 



e^ = <l+2Eh 
Notice the following simple ratios (valid at IPN order) 



1/2 



et E 

- = l+(8-3y)-, 
er c^ 

e, E 

- = l+(8-2y)-, 

er . E 

- = 1 + V— . 



(2.14a) 
(2.14b) 



(2.15) 



(2.16a) 
(2.16b) 
(2.16c) 



(2.17a) 
(2.17b) 

(2.17c) 



In the following paper 112911 we shall need and use the explicit solution of the generalized QK 
binary motion up to 3PN order. 



ni. FOURIER DECOMPOSITION OF THE BINARY'S MULTIPOLE MOMENTS 

A. Peters & Mathews derivation of the Newtonian energy flux 

The method we shall use in this paper is exemplified by the computation of the averaged en- 
ergy flux of compact binaries at Newtonian order using a Fourier decomposition of the Keplerian 
motion yj]. The GW energy flux, say 



r 



^(fr^(/ 



dO. 



d£ 
dtdQ 



GW 



(3.1) 



^ Thus it is sometimes useful to define k' - k/3 which reduces to l/{c^h^) at IPN order. 



where £ is the energy carried in the gravitational waves, reduces at Newtonian order to the standard 
Einstein quadrupole formula^ 

r^' = lfir(t)'iSti (3.2) 

where (N) means the Newtonian limit, the superscript (n) refers to differentiation w.r.t time n times, 
and !■■ is the symmetric-trace-free (STF) quadrupole moment at Newtonian order given by 

/f^=//x<V>. (3.3) 

Here x' is the binary's orbital separation, and the angular brackets around indices indicate the STF 
projection: x'^'x^^ = x'x^ - \^6'h^. Peters & Mathews yj] obtained the expression of the (averaged) 
Newtonian flux for compact binaries on eccentric orbits by two methods. The first method was to 
take directly the average in time of Eq. (13.21 ) using the expression (13.31) computed for the Keplerian 
ellipse; the second method was to decompose the components of the quadrupole moment into 
discrete Fourier series using the known Fourier decomposition of the Keplerian motion (the two 
methods, as expected, agreed on the result). 

In the second method the quadrupole moment, which is a periodic function of time at Newto- 
nian order, is thus decomposed into the Fourier series 



4r^(0 = y Jfe'''^ (3.4a) 

p=~co 

with Jf = r^7jN),->.^, (3.4b) 

ip)^ Jo 2;r 'J 



where i is the mean anomaly of the binary motion, Eq. (12.41) . Since /. . is real the Fourier discrete 



coefficients satisfy (p)lfl^ = (-p)lf^^* {* denotes the complex conjugate). Inserting Eqs. (|3.4|) 
into (|3.2|) we obtain 

-. +CXJ +00 

^(N) ^ i y y (ip ^)3(i „)3 ^ (N) ^ (N) ^,p.,),_ (3_5^ 

5 ^-^ ^-^ (pf (qV 

p=-co q=~-co 

Next we perform an average over one period P which means the average over £ = n{t - tp) which 
is easily performed with the formula 

{e'''')^£^e'''' = Sp^o. (3.6) 

This immediately yields the averaged energy flux in the form of the Fourier series 

^ +00 

(r^^^) = ±yipnf\I^;;)\ (3.7) 



From now on we set c = 1 and G - 1. 



Using dimensional analysis (and the known circular orbit limit) this flux is necessarily of the form 

(^(N)^ = |/(^)V(e), (3.8) 

where v = fj./m and a is the semi-major axis of the Newtonian orbit, and the function f{e) is 
a dimensionless function depending only on the binary's eccentricity e. The coefficient in front 
of (13.81 ) is chosen in such a way that f(e) reduces to one for circular orbits, i.e. when e = 0. Thus 
we have 

The Fourier coefficients of the quadrupole moment are explicitly given by Eqs. (I A3 1) in the Ap- 
pendix |A]below. Remarkably this function admits an algebraically closed-form expression, crucial 
for the timing of the binary pulsar PSR 1913-1-16 HOA . and given by 



1 + 23g2 ^ 37^4 

/(^) = ^^.,2ylP. ■ (3.10) 

The function f(e) is the Peters & Mathews [1] "enhancement" function, so designated because 
in the case of the binary pulsar, which has eccentricity e = 0.617 ■ • ■ , it enhances the effect of 
the orbital P by a factor ~ 11.843. The proof that the series (|3.9I) can be summed up to yield 
the closed-form expression (13.101) is given in the Appendix of 111]. Of course Eq. (13.101) is in full 
agreement with the direct computation of the average performed in the time domain 111], i.e. 

f(e) = :^^^(I^'fH (3-11) 

The method of decomposing the Newtonian moment of compact binaries as discrete Fourier 



series was used in Ref. ||2(J] to compute the tail at the dominant 1 .5PN order. To extend this result 
we need to be more systematic about the Fourier decomposition of the (not necessarily Newtonian) 
source multipole moments. 

B. General structure of the Fourier decomposition 

The two sets of source-type multipole moments of the compact binary system are denoted by 



/z,(0 and J lit) following Ref. 11411] . Here the multi-index notation means L = zi/2 • • • ii, where / is 
the number of indices or multipolarity (which is not to be confused with the mean anomaly {). In 
this Section we investigate the structure of the mass and current moments It and say 7^_i (where 
L- 1 = ?i?'2 • • • ?/-i is chosen in the current moment for convenience rather than L), at any PN order 
and for a compact binary system moving on a general non-circular orbit^. Their general structure 



However the intrinsic spins of the compact objects are neglected, so the motion takes place in a fixed orbital plane. 



can be written as 

hit) = _^:r,[r,r,i;2]x<''-'* {;'■*-'■■■'■'>, (3.12a) 

k=0 

1-2 

J^_i(t) = ^ ^^[r, r, v^] x^''-''v''*'-''-'s''-""''':^v\ (3. 12b) 

k=0 

where e'"^ is the Levi-Civita symbol (such that s^^^ = 1), where x' = y\ -1/2 and v' = dx' /dt = v\ -v\ 
denote the relative position and ordinary velocity of the two bodies (in a harmonic coordinate 
system). In (13.121 ) we pose for instance jc" ■'* = x" • • ■ x'*, and the angular brackets surrounding 
indices refer to the usual symmetric-trace-free (STF) projection with respect to those indices. 

Using polar coordinates r, in the orbital plane (as in Sec. Ill Al) . the above introduced coeffi- 
cients Tk and Qk depend on the masses and on r, r and \?- = f- -^ r^cp^. For quasi-elliptic motion 
we can explicitly factorize out the dependence on the orbital phase (p by inserting x = r cos (f>, 
i/ = r sin 0, and Vx = r cos cp - r^s'mcp, Vy = rs'mcp + r^ cos 0. Furthermore, using the explicit 
solution of the motion (Sec. IIIBI) we can express r, r and v^, and hence the !7^'s and Qk^, as peri- 
odic functions of the mean anomaly £ = n(t - fy), where n = In I P. We then find that the above 
general structure of the multipole moments can be expressed in terms of the phase angle (p, as the 
following finite sum over some "magnetic-type" index m ranging from -/ to +1, 



hit) = y J?lL(^)e""^ (3.13a) 

JL-i{t) = y SL-l(^)e""^ (3.13b) 



m=-l 



involving some coefficients (m)^L and {m)Si-\ depending on the mean anomaly i and which are 
complex (6 C). (Some of these coefficients could be vanishing in particular cases.) The point for 
our purpose is that these coefficients mq periodic functions of £ with period In. As we can see, the 
structure of the mass and current moments Ii and 7l_i is basically the same, but their coefficients 
{m)^L and (m)SL-i will have a different parity, because of the Levi-Civita symbol entering the 
current moment Jl-\- 

To proceed further, let us exploit the doubly periodic nature of the dynamics in the two variables 
A = Ki and £ (as reviewed in Sec. IIIAI) . The phase is given in full generality by Eq. (12.81 ) where 
we recall that W{€) is periodic in i. In the following it will be more convenient to single out in 
the expression of the phase the purely relativistic precession of the periastron, namely A- £ = k£ 
where k = K-l. Inserting the expression of the phase variable into Eqs. (13.131 ) yields many factors 
which do modify the coefficients of (13.131) . but in such a way that they remain periodic in £. Hence 
we can write 



hit) = y I Me''"'', (3.14a) 

7i_i(0 = y Ji-me""'', (3.14b) 



, (in) 
m=-l 



10 



where the coefficients (,n)lLifi) and (m)J^L~i(.^) are 2;r-periodic. Finally, this makes it possible to use 
a discrete Fourier series expansion in the interval £ 6 [0, In] for each of these coefficients, namely 



Idi) = y I Le'"', (3.15a) 

(m) ^-~' (p,m) 

p=-oa 

JL~m = Yj ^^-1^''''' (3.15b) 



(m) ^"^ (p,w) 



with inverse relations given by 

I L = ^ Idi)e-'P', (3.16a) 

(P,m) Jo In (m) 

^2?!- J/ 

J^L-1 = ,^JL-me~'P'. (3.16b) 

This leads then to the following final decompositions of the multipole moments, 

-1 

(P,m) 



hit) = y y J ^,HP^mk)i^ (3_i7^) 

ptt ,t^, (/>■'") 

jL-iit) = y y jl^.c'^p^"''^'. (3.17b) 

^—^ ^—^, (p,m) 
p=-oo m=-\ ' 

Obviously, since the moments //, and 7^-1 are real, their Fourier coefficients must satisfy (p,m)^L = 

(-p-m)J^ I and {p^m)UL-Y = {-p-m)Ui^\- 

The previous decompositions were general, but it is still useful to introduce a special notation 
for the particular case of the Newtonian (N) order, for which the relativistic precession k tends to 
zero. In this case we recover the usual periodic Fourier decomposition of the moments [generaliz- 
ing Eqs. (13.41 )1. with only one Fourier summation over the index p, so that 

+cx> 

lf\t) = y jf^e'''^ (3.18a) 

^^ (p) 

p=-oo 

Jf\{t) = y J,^_7e'^^ (3.18b) 

(p) 

P= — OQ ^^ ' 

The Newtonian Fourier coefficients are equal to the sums over m of the doubly-periodic Fourier 
coefficients in Eqs. (13.171) when taken in the Newtonian limit, namely 



ir = S I T\ (3-19a) 

ip) ,frZi (P''") 



vi-/ m=-l ' 
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IV. TAIL CONTRIBUTIONS IN THE FLUX OF COMPACT BINARIES 

The technique of the previous Section is applied to the computation of the tail integrals in the 
energy flux of compact binaries. Although the computations are efi'ectively done up to the 3PN 
level, the method we propose could in principle be implemented at any PN order. 

A. Expression of the tail integrals in the 3PN energy flux 



As reviewed in the Introduction, the first hereditary term in the energy flux T occurs at the 
1.5PN order and is due to GW tails caused by interaction between the mass quadrupole moment 
and the total ADM mass. At the 3PN order, three kinds of hereditary terms appear: (1) The 
tails caused by quadratic non-linear interaction between higher-order multipole moments with the 
mass; (2) the "tails of tails" due to the cubic non-linear interaction between the tail itself and the 
mass; (3) a particular "tail-squared" term arising from self-interaction of the tail'°. 

In the equations to follow, we list the expressions for all these hereditary tail terms. They are 
given as non-local integrals over the source multipole moments of the system Iij{t), ^ijkilli ■■■ and 
Jijit), ..., where we use the specific definition of the PN source moments given in Ref. [14 ill . Thus 
the energy flux T^ defined by Eq. (13.11 ) can be split at 3PN order into 



r 



■(3PN) _ 



— -* inst + Jh 



heied > 



(4.1) 



where the "instantaneous" part, which depends on the source moments at the same instant (say t), 
reduces at the Newtonian order to the Einstein quadrupole moment flux T^^^ given by Eq. (13.21 ). 
On the other hand, the "hereditary" part reads 



-* hered — -'tail + -* tail(tail) + -* (tail)^' 



where the quadratic-order tail integrals are explicitly given by (see Ref. [13 ill 'j^' 



r. 



tail 



f 

Jo 



189 '^* 



ir^it) 



+ 



64M ,„ 



r(3)/ 



f 

Jo 



dTllf{t-T) 



dri'Sit-T) 



drjf/it-T) 






(4.2) 



(4.3) 



'" Recall that the hereditary character of the non-linear memory integral I 111 lid 
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H 



is that of a time anti- 
derivative in the waveform (i.e. the radiative moments). Thus the non-linear memory becomes instantaneous in 
the energy flux, which is made out of time derivatives of the radiative moments, and will be included into the 
instantaneous terms computed in I29i1 . 
" For convenience we do not indicate the neglected PN terms, e.g. 0(c^"). All equations are valid through the aimed 
3PN precision. In the companion paper 12911 we shall restore all powers of 1/c (and G). 
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while the cubic-order tails (proportional to M^) are 

4M2 



r. 



tail(tail) 



5 '?« 



Jo 



drifXt-T) 



57 



In^l — l + ^ln|^| + 

2n^ 



70 \2r, 



124627 
44100 



J(<im\f - 



AM^ 



Jo 



dTl^{t-T) 



'Ai^yi 



(4.4a) 
(4.4b) 



In these expressions recall that M is the conserved mass monopole or total ADM mass of the 
source. The first term in (|4.3I) is the dominant tail at order 1.5PN while the second and third 
represent the sub-dominant tails both appearing at order 2.5PN. The higher-order tails are not 
given since they are at least at 3.5PN order (see |12fl for their expressions). The two cubic-order 
tails given in Eqs. (14.41) are both at 3PN order. 

The constant tq scaling the logarithms in the above tail integrals has been defined to match 
with the choice made in the computation of tails-of-tails in Ref. lll2n . This is the length scale 
appearing within the regularization factor (r/ro)^ used in the multipolar moment formalism valid 
for general sources 114 ill . Note that ro is a freely specifiable constant entering the relation between 
the retarded time in radiative coordinates [used in Eqs. (I4.3I) - (I4.4I) 1 and the corresponding time in 
harmonic coordinates. Hence vq merely relates the origins of time in the two coordinate systems 
and is unobservable. 

We shall compute all the tail and tail-of-tail terms (|4.3I )-( |4.4|) {i.e. up to the 3PN order] averaged 
over the mean anomaly I. Together with the instantaneous terms reported in the next paper [l29j] 
we shall obtain the complete expression of the 3PN energy flux. It is clear from Eqs. (|4.3I )- (|4.4|) 
that all the terms necessitate an evaluation at the relative Newtonian order except the mass-type 
quadrupolar tail term - first term in (14.31 ) - which must crucially include the IPN corrections. We 
start with all the terms required at relative Newtonian order and then tackle the more difficult IPN 
quadrupolar tail term. 



B. Tails at relative Newtonian order 



As a warm up, we consider the mass-type quadrupolar tail term in the energy flux, the first term 
in Eq. (|4.3I) . but given simply at the relative Newtonian order, namely^^ 



^ mass quad't^" 



( 



AM (3).^. r 



(5) 
dT h^ 



(f-T) 



^"fe)"n 



), 



(4.5) 



where the brackets () refer to the average over the mean anomaly i as defined by Eg. (|3.6I) . The 
term (14.51) was already computed using a Fourier series at Newtonian order in Ref. 1201]; note that 
the method of Urn is valid only for periodic motion and thus is applicable only at the Newtonian 
level. In this Section we recover the Newtonian result of ||20|] . 

The Fourier decomposition of the Newtonian quadrupole moment was already given in general 
form by Eqs. (13.41) . We insert that decomposition into the flux (14.51) and we evaluate the tail 
integral by using the fact that if i{t) = n(t - tp) corresponds to the current time t, then clearly 
£(t - r) = £(t) - nr corresponds to the retarded time t - r. Next we perform the average over the 



'^ We shall compute this term at IPN relative order in Sec. lIVDl 
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current value £{t) with the help of the formula (13.61) . The result is 



(rl 



(N) 
mass quad 



)taii = -^ ) Apnfll 



f 

Jo 



+ 00 



^"'2^)"T^ 



(4.6) 



It remains to handle the last factor in (14.61) which is the tail integral in the Fourier domain, and 
which is computed using the closed-form formula 



dre'^'ln —\ = 



2ro 



(t12 



^signicr) + i(ln{2\cr\ro) + C) 



(4.7) 



where cr = pn, sign(cr) = ±1 and C = 0.577 • • • denotes the Euler constant. Inserting Eq. (I4.7|) 
into (|4.6I) . we check that the imaginary parts cancel out, and the result reduces to 



^ mass quad 



4nM^ , 
)taii = ^y(pn)'lJ,. 



(N)|2 



(4.8) 



p=i 



Observe that the range of p's corresponds to positive frequencies only. Eq. (14.81) agrees with the 
result of [|2(3n and can interestingly be compared with the expression of the Newtonian part of 
the averaged flux (quadrupole formula) as given by Eq. (13.71) . Although Eq. (14.81) is expressed 
in terms of the relatively simple Fourier series (14.81) [unlike for the case of the IPN quadrupole 



tail in Sec. IIVDI which will turn out to be substantially more intricate], it has to be left in this 
form since no analytic closed-form expression can be found for the infinite sum of these Fourier 
components 12(311 . This is in contrast with the quadrupolar Newtonian flux (13.71) which does admit a 
closed-form expression [recall Eq. (13.101) 1. In Sec.|V]we shall further proceed following Ref. [|20D 
by expressing Eq. (14.81) in terms of a new "enhancement" factor depending on the eccentricity and 
which will be computed numerically. 
Let us stress that the result (14.8 



and all similar results derived below are "exact" only in a 

PN sense. Indeed we have formally replaced inside the tail integral the orbit of the binary at any 
earlier time t - r hy its orbit at the current time t, thereby neglecting the effect of the binary's 
adiabatic evolution by radiation reaction in the past. As a result there should be a remainder term 
in (14.81 ), given by the order of magnitude of the adiabatic parameter ^lad = d>/oj^ associated with 
the binary's inspiral by radiation reaction. Indeed, we know lIllL 12011 that the replacement of the 
current motion inside the tail integral is valid only modulo some remainder O (^rad) or, rather, 
O (^lad In^rad)- In tcrms of a PN expansion such remainder brings a correction of relative 2.5PN 
order which is always negligible here (indeed the tails are themselves at 1.5PN order so the total 
error due the neglect of the influence of the past in the tails is 4PN). 



The other tail integrals, given by the second and third terms in Eq. (14.31) . are evaluated in exactly 
the same way. With the PN accuracy of the present calculation these integrals are truly Newtonian 
so the mass octupole moment /,y^. and current quadrupole moment Jtj are required at Newtonian 
order only. For simplicity, we do not add a superscript (N) to indicate this because there can be no 
confusion with other results. We thus need to evaluate the time-averaged fluxes 



\J mass Oct /tail ~ ^ ^ or\ iik^ ' 



64M 



r(3)/ 



\-'cun- quad /tail — \ <^ Jij\J) 



drlflit-T) 



^"'2^)^1 



), 



f 

Jo 



drJlfit-T) 



'"'^)^s 



(4.9a) 
(4.9b) 
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Inserting the Fourier decomposition of the moments, performing the average using Eq. (13.61) and 
using the integration formula (14.71) immediately results in 



(Tmass octXaii = ^ Y ip nf \ lijut (4. 10a) 

189 ^ (p) 

64;rM ■^ ^ , 

(^.u. quad)tan = —T^ ) ip «)' \Jij?. (4. 10b) 

In Sec. |V]we shall have to provide some numerical plots for the eccentricity-dependent enhance- 
ment factors associated with Eqs. (14.101) . since they cannot be computed analytically. 



C. Tails-of-tails and tails squared 

We have seen that at the 3PN order {i.e. 1.5PN beyond the dominant tail) the first cubic non- 
linear interaction, between the quadrupole moment /,y and two mass monopole factors M, appears. 
Following Eqs. (14.41) we thus have to compute the "tail-of-tail" contribution. 



and the so-called "tail squared" one. 



-I- 



124627 
44100 



), (4.11) 



2roj 12 



(n,nf) = (-^[l drl]f(t-T) 



). (4.12) 



Both contributions are evaluated at relative Newtonian order, inserting the Fourier decomposition 



of the Newtonian quadrupole moment (13.41) [suppressing the notation (N) for simplicity]. The new 
feature with respect to the previous computation is the occurrence of a logarithm squared in the 
tail-of-tail integral (14.111) . The integration formula required to deal with this term is [compare with 
Eq.diJ])] 



Jo 



^^dre-nn^i^] = ^l^-- 



y2roj cr [6 
and with this formula, together with (14.71) . we obtain the result 



^sign((r) + i(ln(2|(r|ro) + c) 



(4.13) 



4M2-^ „ i^),(n'^ , x2 57/ ^ 1246271 

(^Tta^KtaU)) = -^YiP ^) \I^I^^\ W - 2[ln{2p n r,) + c) + ^(ln(2;? n r,) + c) - -^^^\ ■ 

(4.14) 
On the other hand the tail squared term is readily computed with (14.71) and found to be 

(^aii)^) = -^Z^P nf \ip jy + 2(ln(2;, . ,„) + c - jl) | . (4.15) 
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Summing up the two results (14.141) and (14.151) we finally obtain 

,_ 4M^^ o (N^,f2;r2 214 214 1167611 

<?;.«.„,««.> = —Bp ")* li™!' |- - To5 l"<2p« n,) --c-^. (4.16) 

As we can see the contribution from logarithms squared has cancelled out between the two 
terms (I4.14I) - (I4.15I) . Such cancellation is in fact known to occur for general sources UM- We 
observe also that the result ( 14.161 ) still depends on the arbitrary length scale tq. It will be important 
to trace out the fate of this constant and check that the complete energy flux we obtain at the end 
(including all the instantaneous contributions computed in ll29ll ) is independent of tq. 



D. The mass quadrupole tail at IPN order 

Let us now tackle the computation of the mass quadrupole tail at the relative IPN order, namely 

<!r„,assquad)tail = (— llfiO J^ dj ff^ {t - T)[ln| — 1 + -|). (4.17) 

At the IPN order (and similarly at any higher PN orders), we must take care of the doubly-periodic 
structure of the solution of the motion [Sec. Ill All , and decompose the multipole moments accord- 
ing to the general formulas (13.171) . So the IPN mass quadrupole moment /;y entering Eq. (14.171) is 

decomposed as 

+00 2 

^0-(0= y y I ije'^P^"^'^', (4.18) 

p=~co m=-2 

with doubly-indexed Fourier coefficients (p,m)Iij which are valid through order IPN. We can be 
more precise and notice that the harmonics for which m = ±1 are zero at the IPN order, so that 

^■W= y I ^ ,7^'^'"''^'+ Iije""+ lije'^P^^''^'], (4.19) 

but in the following it is more convenient to work with the general decomposition (14.181) . keeping 
in mind that the terms with m = ±1 are absent. As before we insert (|4.18l) into (|4.17l) to obtain 
[after neglecting 2.5PN radiation reaction terms (9 (^rad)] 

(^Tmass quadXail = —^ Y^i H^ip + mkf(p' + m'kf I fj J ij 



^—^ , ip,fn) ip' ,m' ) 

p,p'\m,m' v/-' / v/ 



X (ei(P+P'+('"+'"'»f) r dje--^ip'^'n'k)nryi±\ 



.nl^HJi 



(4.20) 



where the summations range from -00 to +00 for p and p' , and from -2 to 2 for m and m' . 
Evidently the factors {p + mk)^ and (p' + m'k)^ come from the time-derivatives of the quadrupole 
moment. We have explicitly left the last two factors in (14.201 ) as they are, namely the average over 
■£ of an elementary "doubly-periodic" complex exponential, and the Fourier transform of the tail 
integral. 

The expression (14.201) is to be worked out at the IPN order. Since the relativistic advance of 
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the periastron k is already a small IPN quantity, the first thing to do is to evaluate (14.201) at linear 
order in k {i.e., neglecting 0{k^) which is at least 2PN]. Afterwards we shall insert the explicit 
expressions for the IPN quadrupole moment and ADM mass. We provide here the necessary 
formulas for performing the linear-order expansion in k of the last two factors in (14.201) . The 
average we perform is over the orbital period (time to return to the periastron) and so is defined by 

/^i(p+mk)e\ — I ^^i(p+mk)t (4 21) 



Jo 2;r 



Using the fact that m fc «c 1 since we are in the limit where k ^> (hence p + mk is never an 
integer unless k = 0), we readily find 

{m \ 

-k if ;? ^ I 
P \ + 0{k^). {A. 22) 

I + inmk ifp = OJ 

This result depends only on whether p is zero or not, and is true for any integer m, except that 
when m = the result (14.221 ) becomes "exact" as there is no remainder term 0(k^) in this case. 

On the other hand, to compute the tail integral given by the last factor in Eq. (14.201) . we expand 
it at first order in k, obtaining thereby 

P jTe'('-'"^)"Mn(-^) = (l - — ) P Jre'^'"Mn(-^]-i^ +0(^2), (4.23) 

Jo \2ro/ \ p /Jo \2rol p^n 



and we apply for the remaining integral in (|4.23|) the formula (|4.7I) . 

With Eqs. ( 14.221) and (14.231 ) in hand we can explicitly work out the tail expression ( 14.201 ) at first 
order in k (the extension to higher order in k would in principle be straightforward). The result 
will be left in the form of the multiple Fourier series (14.201) . into which the results (|4.22|) - (|4.23l) 
have been inserted (we do not try to give a more explicit form for this result which is given by a 
complicated Mathematica expression). In the next Section we shall re-express this series in terms 
of some elementary enhancement functions which will finally be evaluated numerically. 

V. NUMERICAL CALCULATION OF THE TAIL INTEGRALS 

A. Definition of the eccentricity enhancement factors 

We define here some functions of the eccentricity by certain Fourier series of the components 
of the Newtonian multipole moments /[ and 7^ j for a Keplerian ellipse with eccentricity e, semi- 
major axis a, frequency n = 2nlP (such that Kepler's law n^a^ = m holds at Newtonian order). In 
the frame of the center of mass we have f^^ = ^si{v)x^^^ and j''^\ = iJ.si(v)x'^^'^£''-'^"'^x"v'^ where 
yu = mim2/m = vm. Here we pose sj(v) = X'^^ + (-)'X|"^ , where X^ = ^ = Ul + Vl -4v) , 

and -^2 = ^7 = ^1" Vl - 4vj . Let us rescale the latter Newtonian moments in order to make 
them dimensionless by posing 

/f^ = iua's,(v)lL, (5.1a) 

7^ = iua'nsi(v)JL^,. (5.1b) 
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Our first "enhancement" function is of course the Peters & Mathews Uj] function which we have 
already expressed in Eq. ( 13.91) as a Fourier series [and which turns out to admit the analytically 
closed form (13.101) 1. In terms of the Fourier components of the rescaled quadrupole moment Uj 
this series reads 

-. +CX) 

and is such that the averaged energy flux of compact binaries at the Newtonian order reads 

<5-^^^> = yv2xV(e), (5.3) 



where we have defined for future convenience the frequency-related PN parameter x = {mco) 



2/3 



where uj is the binary's orbital frequency defined for general orbits by Eq. (12.71) . Note that in 
Eq. (15.31) which is Newtonian we can approximate cohy n (hence x reduces to m/a). 

Next, we define several other "enhancement" functions of the eccentricity which will permit to 
usefully parametrize the tail terms at Newtonian order. First we pose 

^W^SE"'!!/- (5.4) 

p=\ 

Like for f{e) this function is defined in such a way that it tends to one in the circular orbit limit, 
when e ^ 0. However, unlike for /(e), it does not admit a closed-form expression, and will have 
to be left in the form of a Fourier series. The function ip{e) parametrizes the mass quadrupole tail 



at Newtonian order, in the sense that we have, from Eq. (14.81) . 

(^ri,uJ = fv^-l4..^/V(e)]. (5.5) 

For circular orbits, ^(0) = 1 and we recognize the coefficient An of the 1.5PN tail term (oc jc^''^) 



as computed numerically in Ref. [|42D and analytically in Refs. [|20l . 14311 . The function ip{e) has 
already been computed numerically from its Fourier series (15.41) in Ref. 1201. Here we show the 
plot ofip{e) in Fig.[i|(see Sec. lVBI for details on the numerical computation)'^. 

We next proceed similarly for the 2.5PN mass octupole and current quadrupole tails. We pose 

20 ^"^ 
Pie) = 7TrT7^y/|J/,iP, (5.6a) 

49209^ ip)^ 
p=\ 

+0O 



y{e) = AYp'\Jijt (5.6b) 

Again these functions tend to one when e ^ (as will be checked later) and most probably do not 



'■^ Note that our notation is different from the one in f2Cfl; the function i^Bs(e) there is related to our definition by 
V'Bs(e) - if{e)/f(e). In the present work it is better not to divide the various functions by the Peters & Mathews 
function /(e) entering the Newtonian approximation. 
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FIG. 1 : Variation of (p{e) with the eccentricity e. The function (p{e) agrees with the numerical calculation of 
Ref. 12011 modulo a trivial rescaling with f{e). The inset graph is a zoom of the function (which looks like a 
straight horizontal line in the main graph) at a smaller scale. The dots represent the numerical computation 
and the solid line is a fit to the numerical points. In the circular orbit limit we have (/j(0) = 1. 



admit any closed-form expressions. With their help these tail terms (oc x^^^) of Eqs. (14.91) read 



_ 32 2 5 
\-' mass oct/tail — ^ ^ -^ 

'j'curr quad /tail ~ ~c~ 



n(l-4v)x^'^/3(e) 



16403 
2016 



(5.7) 
(5.8) 



The numerical graphs of the functions jS(e) and y(e) are shown in Fig.[2l 

Two further enhancement factors are then introduced to parametrize the tail-of-tail and tail 
squared integrals (which are Newtonian with the present approximation). The first of these func- 
tions looks very much like the Peters & Mathews function f{e), Eq. (15.21 ). in the sense that its 
Fourier series involves even powers of the modes p. Namely we define 



-. +00 



(5.9) 



Thanks to this even power oc p^ vve find that F{e) can also be computed as an average performed 
in the time domain similar to the one of Eq. (13.111) for /(e). Namely we easily verify that 



F{e) = 



1 



'f(4)f(4)^ 



128 n«^ '^ '^ ^' 



(5.10) 
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FIG. 2: Variation of jS(e) (left panel) and y{e) (right panel) with the eccentricity e. In the circular orbit limit 
wehave/3(0) ^7(0) = 1. 



which can straightforwardly be computed in the time domain with the result that F{e) admits like 
for /(e) an analytic closed form which is readily obtained as 



F{e) = 



1 + 85^2 , 5171 4 , ITM 6 , ^21 8 
6 192 "^ 192 "^ 1024*^ 



(1 - e2)13/2 



(5.11) 



On the other hand we shall need to introduce a function whose Fourier transform differs from the 
one of F{e) by the presence of the logarithm of modes, namely 






p=\ 



(p) 



(5.12) 



One can be convinced that very likely ;^^(e) does not admit any analytic form [hence we name it 
using the Greek alphabet - in contrast to /(e) and F(e)]. Note that;t'(e) has been exceptionally 
defined in such a way that it vanishes when e ^ 0. This is easily checked since in the circular 



(N) 



orbit limit (and at Newtonian order) the quadrupole moment l\. possesses only one harmonic 
corresponding to p = 2 which due to the log term reduces ;t'(e) to zero in this case. In Fig. [3] we 
show the numerical plot of the function ;^^(e) [and also the one for F(e)]. 

With those definitions we find that the sum of tail-of-tail and tail squared contributions obtained 
in Eq. (|4.16|) reads 



kT, 



tail(tail)+(tail) 



2) 



32 



v'x^ 



116761 16 , 1712^ 1712, ,, 

+ — Ti C In (4 u) rn) 

3675 3 105 105 ^ " 



1712 
FW-— ^W 



(5.13) 
The circular-orbit limit can be immediatel y re ad off from this expression and seen to agree with 
Eq. (5.9) in Ref. HI or Eq. (12.7) in Ref. mil. 

Finally we provide the result in the case of the mass quadrupole tail at IPN order. We have seen 
in Sec. IIVE)] that the calculation in this case is much more involved, as the Fourier series (14.201) 
contains several summations, and depend on the intermediate results (14.221) and (14.231) . In addition 
the computation must take into account the IPN relativistic correction in the mass quadrupole 
moment and ADM mass; these are provided in Eqs. (15.161) and (15.171) below. We find that probably 
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FIG. 3: Variation of ;^'(e) (left panel) and F{e) (right panel) with the eccentricity e. In the right panel, 
the exact expression of F{e) given by Eq. (I5.11I) is used. In the circular orbit limit we have ;\'(0) - and 
f (0) = 1. 



there is no simple way {i.e. no simple-looking Fourier series like for instance (15.121) 1 for expressing 
the new enhancement functions of eccentricity which appear at the IPN order. However one can 
check beforehand that the IPN term is a linear function of the symmetric mass ratio v, hence we 
must introduce two enhancement functions, denoted below a and 6. As before we normalize these 
functions so that a{0) = 1 and 9(0) = 1. We have [extending Eq. (15.51) at the IPN order] 



(r. 



mass quad /tail 



)t, 



^y^x'"44nipie,) + nx 



428 , , 178 ^, ■ 
a(et) + -—ve(et) 



21 



21 



(5.14) 



This equation provides the definition of the two enhancement functions a and 6, and we resort to 
the Mathematica computation to obtain them as complicated Fourier decompositions, which will 
then be directly computed numerically using the method outlined in Sec. IV B[ Notice that since 
we are at the IPN level we must use a specific definition for the eccentricity, and we adopted here 
the "time" eccentricity e, entering the Kepler equation (12.1 Obi) in Sec. IIIBI At the IPN order the 
other eccentricities are related to it by Eqs. (I2.17D . On the other hand, the frequency-related PN 
parameter, given by 

x = (mcof\ (5.15) 

crucially includes the IPN relativistic correction coming from the periastron advance K = I + k, 
through the definition o) = nK of Sec. IIIAI All the IPN corrections arising from the formu- 
las (14.221 ) and (14.231 ). the multipole moments M and Ijj, the use of the time eccentricity e, and the 
specific PN variable x, are incorporated in a Mathematica program dealing with the decomposi- 
tion (14.201 ) and used to obtain (15.141 ). The behaviour of the enhancement functions a{e) and 6(e) 
are given in Fig. HI 



B. Numerical evaluation of the Fourier coefficients 

We now describe the numerical implementation of the procedure for the computation of the 
Fourier coefficients of the multipole moments that lead to the numerical plots of the previous 
Section. We focus the discussion on the computation of the crucial coefficients (p,m)^ij at IPN 
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FIG. 4: Variation of a{e) (left panel) and 6{e) (right panel) with the eccentricity e. In the circular orbit limit 
wehavea(O) = 9(0) = 1. 

order which are the more difficult to obtain. The mass quadrupole moment with IPN accuracy is 
given by [compare with the general structure (I3.12al) 1 



/,-,■ = // 1 + 



29 29 



m 



V + — — + - V 

42 14 / r \ 7 7 



8 



x^'xJ^ 



11 11 \ 2 a /) / 4 12 , ,,. ., 



(5.16) 



where x' and v' = dx'/dt are the relative position and velocity in harmonic coordinates, and r = \x'\ 
(like in Sec. IIIBI) . Equation (15.161) is valid for non- spinning compact binaries on an arbitrary 
quasi-Keplerian orbit in the center-of-mass frame (see e.g. 114411 ). Since we investigate tails with 
IPN relative accuracy we need also the relation of the ADM mass M to the total mass m = mi +m2 
at IPN order, 



M = m 



v^ m 
l+vl--- 



(5.17) 



Using the quasi-Keplerian representation of the motion [Sec. IIIBL the dependence of 7,^ on 
x', v', r, v^ and r can be parametrized in terms of the eccentric anomaly u. However, as explained 
previously we require Iij(-£) in the time domain to proceed. The steps of our numerical implemen- 
tation scheme can be summarised as follows: 

1. To begin with, each component of the IPN mass quadrupole is expressed in terms of 
the quasi-Keplerian parameters using Eqs. (I2.10I) - (I2.12I) . The components of the mass 
quadrupole are now functions of the eccentric anomaly u, and are parametrized by the mean 
motion n and by one of the eccentricities which is chosen to be e, - the "time" eccentricity 
in Kepler's equation (|2.10b| ). ^'^ 

2. We next invert, numerically, the equation for the mean anomaly ^ = m - e, sin m to obtain the 
function u(£). This can be done either by using the series representation in terms of Bessel 



The semi-major axis Ur and the other eccentricities e,- and e^ are deduced from n and et using Eqs. (I2.14b -( l2.17b . 
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functions, 



+00 ^ 
= £ + 2Y^-J,{set)sm{sl\ 



(5.18) 



s=\ 



or numerically by finding the root of £ = u - e, siuM. The latter is a more efficient and 
more accurate method and we employed it in this work (we used the FindRoot routine in 
Mathematica). In this case we generated a table of 20000 points of u and i between 
and In (for each value of e^). The above inversion enables us to re-express all functions of 
the eccentric anomaly u as functions of the mean anomaly L If required, a more accurate 
implementation for solving Kepler's equation along the lines of [45] can be used in the 
future. 



3. One needs to be careful in dealing with the u dependence of V in Eq. (12.121) to avoid the 
discontinuity there. To this end it is best to use 



/ Ba, smM \ 

Viu) = u + 2 arctan 

W - Bj, cosm/ 



(5.19) 



where jS<^ = [1 -(1-ei,) '-]/e^. By this process, we thus have in hand the Fourier coefficients 
{m)Iij{£) defined in Eq. (13. Hal) as explicit (numerical) functions of i. 



4. Recall that these functions also have a dependence on the mass ratio v and the PN parameter 
X defined by (m oSf^^ where to = nK. To avoid assuming numerical values for v and x and 
hence to preserve the full generality of the result, we split the function (,„) J,y into 



Iij(e,et,v,x)= I^(e,et) + x 



(m) 



(m) 



-10/ 



11 



(m) (m) 



ie,ed 



(5.20) 



Notice that we have neglected the terms higher than IPN in writing the above expression. 
Now the various („j)J/'-^ are only functions of i and e,. We evaluate the Fourier coefficients 
of these terms separately in the next step of the procedure. 

5. For a fixed value of e,, we can straightforwardly get the plot of (m)J^ff versus £. Equivalently, 
one can also write the Fourier decomposition of („,) J°°(^) as 



(m) 



^—^ (n n 



P=-° 



m ipe 

ij '^ 

(p,m) 



(5.21) 



Now we seek a numerical fit to Eq. (15.201 ). in powers of e^''^ , to extract out the coefficients 

(p,m)-^ ij 



^I?? . We do the same for different values of e, and for 



(PM)Iij and (p,m)Iij ■ 



6. The fitting procedure mentioned above can be implemented either starting with the STF 
moment /,y or the non-STF projected one. The expressions will be different in these two 
cases as for the first case the zz component of the moment is not equal to zero by definition 
[since I^- = -{hx + Iyy)\ whereas for the latter case the zz component is zero due to planar 
motion. This provides a simple algebraic check on the numerical calculation. 

7. Instead of using the basic multipole moment as the starting function {e.g. /,y), we find that 
using the leading time derivative {i.e. l\f) improves the numerical convergence of the results 
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because one deals with lower derivatives of the basic function. This is very helpful for higher 
values of eccentricity. 

8. Substituting the Fourier coefficients into Eq. (14.201) one can generate the numerical values of 
the averaged energy flux (Tmafu quad) for the different values of e,, and hence get the numerical 
values of the enhancement functions, and most importantly of the IPN ones a(e,) and 6{et) 
defined by (15.141) . The plots of these functions reported in Sec. IV Al readily follow. 

We have just described the procedure for the most difficult IPN quadrupole tail yielding the 
computation of a(e,) and 6{et). This procedure is quite general, and provides a method which 
could be extended to higher FN orders. However, at the Newtonian order it is in fact much more 
efficient to make use of the well-known Fourier decomposition of the Keplerian motion. Using 
this we can derive the components of the multipole moments (at Newtonian order) as a series of 
combinations of Bessel functions. Then it is a very simple matter to compute numerically the 
associated "Newtonian" enhancement functions [namely the functions (fie), /3(e), y{e) and xi.^) 
defined in Sec. IV AH . For the convenience of the reader we give in Appendix |A] all the expressions 
for each of the components of the required Newtonian moments \lf^\ 4T ^^'^ -^T^^ ^^ ^ series of 
Bessel functions. We have used them to compute numerically the enhancement functions ip{e), 
yS(e),r(e)andAr(e)'^ 

VI. THE HEREDITARY CONTRIBUTION TO THE 3PN ENERGY FLUX 

A. Final expression of the tail terms 

Based on the treatment outlined above of a numerical scheme for the computation of the orbital 
average of the hereditary part of the energy flux up to 3PN, we finally provide the complete results 
for the numerical plots of the dimensionless enhancement factors. It is convenient for the final 
presentation to redefine in a minor way the "elementary" enhancement functions of Sec. IV A[ 
which were directly given by simple Fourier decompositions. Let us choose 

13696 16403 112 
Me) = a(e) B(e) y(e), (6.1a) 

1424 16403 16 

ne) = e(e) + B(e) + y(e), (6.1b) 

^^ 4081 ^ 12243^^^ 1749^^ ^ ^ 

59920 
K(e) = F(e) + -——x{e). (6.1c) 

11d/d1 

Considering thus the 1.5PN and 2.5PN terms, composed of tails, and the 3PN terms, composed of 
the tail-of-tail and the tail-squared terms, the total hereditary contribution to the energy flux (14.21) 
when averaged over £ (and normalized to the Newtonian value for circular orbits) finally reads 



(^Thered) = ^ v" x' {AtT x''^ cpie^) + 71 x"^ 



8191 ^ ^ 583 , ■ 

met) vi{e,) 

672 ^^ '^ 24 ^^ '^ 



'^ On the other hand, for the Newtonian tail terms, we could proceed exactly in the same way as for the IPN term, 
following the steps 1-8. We have verified that both methods agree well. 
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FIG. 5: Variation of tfr(e) (left panel) and ^(e) (right panel) with the eccentricity e. The inset graph is a 
zoom of the function (which looks like a straight horizontal line in the main graph) at a smaller scale. The 
dots represent the numerical computation and the solid line a fit to the numerical points. In the circular orbit 
hmit we have (A(0) ^ ^(0) ^ 1. 
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(6.2) 



In this result all the enhancement functions reduce to one in the circular case, when e, = 0, so the 
circular-limit is immediately deduced from inspection of Eq. (16.21) . and is seen to be in complete 
agreement with Refs. [|l2l . bill . The function F(et) is known analytically, and we recall here its 
expression. 



F(e,) = 



\ + §5 2 _L 5171^4 _L 1751^6 , 297 



+ 



192 ^t ^ 192 •=/■ ^ 1024 "^f 



(6.3) 



(1 - e?)i3/2 

However the other enhancement functions ip{et), il/{et), ^(e?) and /c(e,) in Eq. (16.21) (very likely) do 
not admit any analytic closed-form expressions. We have explained in Sec. IV B I the details of the 
numerical calculation of these functions. We now present the numerical plots of the final functions 
i/^(e,), ^(e,) and /c(e;) in Figs. |5]-l6]as functions of the eccentricity e, [recall that the function (p(et) 
has already been given in Fig.[T| '^. 

As seen from Eq. (16.21) the final result depends on the constant ro at the 3PN order. Let us 
understand in bit more detail the occurrence of this constant. We first recall from Ref. [|l2ll that 



the dependence on the constant ro of the radiative quadrupole moment at infinity, say Uij, arises 
precisely at the 3PN order, and comes exclusively from the contribution of tails-of-tails (i.e. the 
cubic multipole interaction M^ x 7^). It is explicitly given by 



Uij(t) 



4?^(0 + 



-I- 



214 
T05 



M^ llfit) In ro + 



(6.4) 



in which we have indicated that Uij simply reduces to the second time derivative of /,y at the 
Newtonian order, and where we show the only term which depends on the constant ro; such a term 
appears at 3PN order and turns out to be proportional to the fourth time derivative of /,y. The dots 
in Eq. (16.41) denote many terms which do not depend on ro. From (|6.4I) it is then trivial to deduce 



' The numerical results used for the figures 1-6 are available in the form of Tables on request from the authors. 
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FIG. 6: Variation of K{e) with the eccentricity e. In the circular orbit hmit we have k{0) = 1. 
that the corresponding dependence on ro of the averaged energy flux at 3PN order must be 



(3PN)\ _ 



/<Tr-(Ji-iN)\ _ 



Iku'^') . 



1 



r(3)7(3)\ 



= 5<^^^0 + 



.|^M^(/-Oln.o. 



(6.5) 



Now we can take advantage of the fact that inside the operation of averaging over £ [denoted by () 
and defined by (14.211) 1 one can freely operate by parts the time derivatives. Hence, we can write 



7(3) 7(5) \ 



7(4) 7(4) \ 



that (/.;./..) = -(/../..) and so we arrive at the result 

U 'J 'J 'J 



(r 



(3PN)\ _ 



7(3) 7(3) \ 






525 ^ 'J 'J ^ " 



(6.6) 



The factor of In ro in Eq. (16.61) looks like a "quadrupole formula" but where the third time derivative 
of the moment would be replaced by the fourth one. Notice that the above expression has been 
computed for general radiative-type moments and is true for any PN source, in particular for a 
binary system moving on an eccentric orbit. Therefore the dependence on Inro found in (16.61) 
should perfectly match with the one we have obtained in Eq. (16.21) . Thus, comparing with (16.21) . 
one readily infers that the function F(e,) in the case of an eccentric binary must necessarily be 
given by the components of the quadrupole moment in the time domain as 



F{e,) 



128 v^x^^ 'J 'J ' 



(6.7) 



This prediction is perfectly in agreement with our finding for the function F{e^ in Eq. (I5.10|) 
(indeed, since we are at leading order, M reduces to m, e^ agrees with e, to equals ri). We have 
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therefore confirmed the correctness of the dependence upon ro of Eg . (I6.2h . 

We already know from the study of the circular-orbit case (cf. 113 ill ) that the dependence on 
ro is cancelled out with a similar term contained in the expression of the source-type quadrupole 
moment Ijj at 3PN order. This cancellation must in fact be true for general sources, and has 
been proved on general grounds in Ref. I112n . It will therefore give an interesting check of our 



calculations when we show in the companion paper Il29l1 that the cancellation of rg occurs for 
general eccentric orbits. 

To finish let us provide here the expressions of our final enhancement functions at the first order 
in e^ when e, -^ 0. These expansions will be useful in the following paper 112911 . when we compare 



the perturbative limit of the complete energy flux at 3PN order (including all instantaneous terms) 
with the result of black-hole perturbations. Note that those expansions are obtained analytically. 
For the functions which are Newtonian we can either use the Fourier coefficients in the AppendixlAl 
and expand them at first order in ej or follow the general procedure explained in Sec. IV Bl for the 
relevant moments but expanding Eq. (15.181) to only first order in e^, namely, 



t ■ 

2 



u = i + etsm£+^sm2e + 0[e^,). (6.8) 

Concerning the two IPN functions [i/^(e,) and ^(e^)], on the other hand, we obtain them directly 
using the latter procedure. We find 

Vfe) = l+^'!?+0(<.;), (6.9a) 

t*fe)= l-^e?+0(ef), (6.9b) 

ffc)=l + ^^.?+0(.?), (6.9c) 

,(„) = 1, (62 ,4613840 24570945 ),. 

^ '^ \3 350283 1868176 ) ' ^ '^ ^ 

and of course [since this is immediately deduced from Eq. (16.31) 1 

F(e,) = l + ^e^+0(et). (6.10) 

We have checked that the numerical results of Figs{Tl [5] and [6] agree well with Eqs. (16.91) in the 
limit of small eccentricities. 



B. Conclusion and future directions 

The far- zone flux of energy contains hereditary contributions that depend on the entire past 
history of the source. Using the GW generation formalism consisting of a multipolar post- 
Minkowskian expansion with matching to a FN source, we have proposed and implemented a 
semi-analytical method to compute the hereditary contributions from the inspiral phase of a bi- 
nary system of compact objects moving on quasi-elliptical orbits up to 3PN order. The method 
explicitly uses the IPN quasi-Keplerian representation of elliptical orbits and exploits the dou- 
bly periodic nature of the motion to average the fluxes over the binary's orbit. Together with the 



instantaneous contributions evaluated in the next paper ll29ll . it provides crucial inputs for the con- 
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struction of ready-to-use templates for binaries moving on eccentric orbits, an interesting class of 
sources for the ground based gravitational wave detectors LIGO/Virgo and especially space based 
detectors like LISA. 

The extension of these methods to compute the hereditary terms in the 3PN angular momentum 
flux and 2PN linear momentum flux is the next step required to proceed towards the above goal 
and is currently under investigation. The extension to compute the 3.5PN terms for elliptical orbits 
is currently not possible due to some as yet uncalculated terms in the generation formalism at this 
order for general orbits. It would also require the use of the 2PN generalised quasi-Keplerian 
representation for some of the leading multipole moments. 
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APPENDIX A: FOURIER COEFFICIENTS OF THE MULTIPOLE MOMENTS 

In this Appendix we provide the expressions of the Fourier coefficients of the Newtonian mul- 
tipole moments in terms of combinations of Bessel functions. We decompose the components of 
the moments as Fourier series, 

+00 

If\t) = y iPe^^^, (Ala) 

^~^ (p) 

p=-ca 

jT-\(t) = E ^L-T^'^^ (Alb) 

(p) 



where the Fourier coefficients can be obtained by evaluating the following integrals 

ir = ^ r^^/fCO^-'^^ (A2a) 

(P) 2n Jo 

Jl^i = ^ rde/^_\(t)e-P^. (A2b) 

(P) 2;r Jo 



For the mass quadrupole moment at Newtonian order we have 



17 



/r = ^g + 5<K.(P^.) 



+ 1-^^? - gen(-^p-i (P(^t) + Jp+i (pe,)) 



Note that the Fourier coefficients we provide are for normalized multipole moments as defined in Eqs ( I5.1al l-( l5.1bb . 
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+ 1 ^ + 4^? 1(^2 (P(^t) + Jp+2 (pet)) 

+ i—^er + — en (^3 (P(^t) + Jp+3 (ped), (A3a) 

IP = -i ^l-e^i^et (-Jp-i (pet) + Jp^i (peS) 
+ 1 -^ - 4^?) (-^^+2 (pet) - Jp-2 (pet)) 
+-^et {jp+3 (pet) - Jp-3 (pe,)) >, (A3b) 



lfy' = \2-^t\Mper) 



+ \o^t+ 4^f I (Vi (P^t) + Jp+i (pet)) 

-J (jp-2 (pet) + Jp+2 (pet)) 

+ 1 g e, - — ^f I ( ^3 (P^t) + Jp+3 (pet)), (A3c) 



r/'-i-l-le^MP.) 
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For the mass octupole moment we find 



+ 1 2^' + g^?) (Vi (P^t) + Jp+i (pet)) 

--e] [jp^2 (pet) + Jp+2 (pet)) 

+—e^ (7p_3 (pet) + Jp+3 (pet)). (A3d) 



ifr =-#'^T<l^''<'-'' 



7-(N) 

(p) 



(3 21 1 1 \ 

"40 " 20"' " 40"') (-^^'-^ ^^''^ ^ -^^^^ ^^''^) 

^ ( 20^' + ^^n (^2 (pet) + Jp+2 (pet)) 

(13 3 \ 

- g - ^^? + ^^n (^3 (P^f) + -^^+3 (pet)) 

+ ij^et - —e^j (jp-4 (pet) + Jp+A (pet)) \, (A4a) 
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+^e^ (-Jp-4 (pet) + Jp^A {pet)) L (A4f) 

Finally, for the current quadrupole moment, 

^,f = -\J]T~ih,eJp{ped 
(P) 4 > ^ 

-- (l + e]) {jp^i {pet) + 7p+i {pet)) 

+-et (jp-2 {pet) + Jp+2 {pet)) \, (A5a) 

J^/f = ^ (l - e?) |(^p+i {pet) - Jp-x {pet)) 

--et (jp+2 {pet) - Jp-2 {pet)) \- (A5b) 
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